How to move a black hole without excision: gauge conditions for the numerical 

evolution of a moving puncture 
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Recent demonstrations of unexcised black holes traversing across computational grids represent 
a significant advance in numerical relativity. Stable and accurate simulations of multiple orbits, 
and their radiated waves, result. This capability is critically undergirded by a careful choice of 
gauge. Here we present analytic considerations which suggest certain gauge choices, and numerically 
demonstrate their efficacy in evolving a single moving puncture black hole. 
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I. INTRODUCTION 

Gravitational waveforms from binary black hole merg- 
ers can only be obtained through 3D numerical relativity 
simulations of the full Einstein equations. Such simula- 
tions have proven to be very challenging. Improvements 
have often come from new formulations or a new choice 
of gauge. For example, a corotating gauge led to the first 
evolution of a binary through one full orbit Q] as well as 
through a plunge, merger, and ringdown 0. 

Recently, the authors of this paper Q and another re- 
search group |4| independently developed the capability 
to numerically evolve freely moving black holes, i.e. black 
holes that are in no way fixed to any position in the com- 
putational grid, without using excision techniques for the 
black hole interiors. These efforts led to the simulation 
of record-breaking numbers of orbits, among other ac- 
complishments . Both of these research programs em- 
ployed gauge conditions that differed not only from those 
prescribed in prior literature, but also from each other. 
A careful choice of gauge is important, particularly for 
moving black holes, as it can determine whether a sim- 
ulation is stable or becomes irrevocably corrupted by a 
runaway build up of numerical error. There are various 
ways in which a choice of gauge can affect the stability 
and accuracy of a numerical simulation; for example, a 
poor choice of gauge might cause large gradients in the 
fields which in turn engender large finite differencing er- 
ror. 

In this paper we investigate a group of gauge choices 
analytically, and numerically demonstrate their effects 
on moving black hole simulations. These gauge choices 
are framed in the context of a conventional, 3+1, con- 
formal formulation of Einstein's equations known as the 
Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formula- 
tion 0, 11,0- After a brief exploration of the shortcom- 
ings of a commonly used gauge choice in ^ we calcu- 
late in section ITlI Al thc characteristic speeds for various 
gauges and then show the numerical behavior of some 
of them in section IIII Bl In section IIVI we discuss the 
properties of the most promising gauge choices and sug- 
gest further improvements. Conclusions are presented in 
section Ivl 



II. ORIGINAL 1+LOG SLICING AND 
GAMMA-DRIVER SHIFT CONDITIONS 



We consider first a gauge introduced in [T( 
non-moving punctures without excision: 

dta = —2aK 
d t [3 l = \a^ n B l 
d t B l = d t r-r/B l 



for evolving 

(1) 
(2) 
(3) 



where a is the lapse, (3 l is the shift, K and f 1 are the 
usual BSSN variables, \&q is the initial conformal factor 
($| £ in H3), B l is an auxiliary variable to make the 
shift equation first order in time, and 77 is a constant 
damping factor. This gauge has proven to be very effec- 
tive at evolving a single non-moving puncture for arbi- 
trarily long durations, and in some circumstances it can 
also be applied to binary punctures. Here the condition 
on the lapse, a variation of the Bona-Masso slicing condi- 
tion [lflllli, is the standard "1+log" singularity-avoiding 
slicing condition. Of particular interest is the "Gamma- 
driver" condition for the shift, which is designed to ul- 
timately "freeze Gamma", i.e. drive dtT to zero as the 
physics of the spacetime also evolve towards quiescence. 

For n > 0, ^g"™ = at the puncture, and thus this spe- 
cific gauge will ensure that f3 l will not evolve there. As 
the motion of the puncture is determined by the shift, the 
puncture will remain motionless in this gauge since the 
shift is initially zero. For black holes that are physically 
dynamical the coordinates may become highly distorted. 
For example, in the case of binary black holes, the separa- 
tion between the black holes as determined by the proper 
distance from horizon to horizon in some spatial slice 
must decrease with time. This implies that components 
of the 3-metric will approach zero between the black 
holes. Then the evolved part of the BSSN conformal 
exponent <j) must either approach —00, or components of 
the conformal 3-metric 7^ must approach zero. In the 
latter case the inverse 3-metric must diverge in order to 
maintain a unit determinant. Thus, in the fixed puncture 
approach, arbitrarily large, ever-increasing fields and gra- 
dients around the black holes seem likely to develop and 
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cause numerical difficulties. Indeed this can be demon- 
strated even in the case of a single moving black hole. 

We will use as a test case a single black hole, given by 
Bowen-York puncture initial data [12j (as computed by 
the elliptic solver AMRMG with unit puncture mass 

and momentum such that its physical speed should be 
half the speed of light. We evolve this data with our 
Hahndol codeQ, which uses the usual, conformal BSSN 
formulation of Einstein's evolution equations on a cell- 
centered numerical grid, with 4 th -order spatial differenc- 
ing and 4 th -order Runge-Kutta time-integration. The 
initial puncture position is at coordinates (— 3M, 0,0) 
(where M is the puncture mass) and the momentum is 
in the positive x-direction. The x-axis is between grid 
points, which are a distanced M/16 apart; all data pre- 
sented here has been interpolated onto the x-axis for plot- 
ting purposes. f3 l and B l are initialized to zero, while the 
lapse is precollapsed with the profile ^q 2 as suggested in 
|l0j| (and recommended in jj| for moving punctures). The 
damping parameter 77 is 2, unless otherwise stated. We 
will assess the performance of each gauge by the evolving 
behavior of the quantities a, f3 x , 4>, and T x ; in particular 
the peak in cj> will roughly indicate the position of the 
puncture and any extreme gradients in T x = —dj^ JX will 
tend to indicate problems with the gauge. 

In this regard it should be noted that T l is completely 
controllable by the gauge condition, and so any unde- 
sirable features observed in T l are in principle avoidable 
via a better choice of gauge. In particular, the specific 
"Gamma-freezing" condition T l = may be desirable, as 
this elliptic generalization of isotropic coordinates pro- 
posed by Dirac is expected to be non-pathological [l4lll5| . 
So, among hyperbolic Gamma-driving conditions, those 
that result in smaller values of T l and its derivatives 
might be preferred. 

In Fig. ^ evolving with the original gauge given in 
Eqs. (lll.'it . we see that T x grows very large and <j> devel- 
ops sharp features. Shortly thereafter the inverse confor- 
mal metric diverges, which is our criterion for stopping 
the run. This failure motivates investigation of moving 
punctures. 

If we wish the punctures to move we must allow P l 
to evolve freely away from zero and thus we must take 
the conformal factor out of the shift equation (setting 
n = Q in Eq. J2J). The results, shown in Fig. |3 are not 
much improved. The puncture is now free to move but 
non-propagating features at x = — 3M destabilize the 
simulation. Clearly there is a problem with zero-speed 
modes. The extended, collapsed region of the lapse points 
to a particular difficulty in this regard, and suggests an 
improvement. 

The factor of a in Eq. (J2J should also be removed, for 
even if it were not precollapsed it would still be driven 
nearly to zero around the puncture by the 1+log slic- 
ing condition, and therefore may retard the evolution of 
f3 l . When this factor is removed, the results are not as 
catastrophic as before, but the simulation is plagued with 
noise. Fig.|3|shows, in particular, a tendency for the min- 
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FIG. 1: Results of the non-moving puncture gauge dta = 

-2aK, dtp = |a*o n B i , d t B i = F - ■qB i at time t = 17M. 
The conformal metric becomes singular by t = ISM. 
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FIG. 2: Results of the gauge dta = -2aK, dtp 1 = \aB\ 

dtB 1 = Y % — r\B % at time t = 17M. Sharp features around 
x = — 3M fail to propagate. The conformal metric becomes 
singular by t — 18M. 



imum in the lapse to lag behind the puncture position. 
The above results suggest that a careful study of the 
propagation speeds in this system might be helpful. 



III. CHARACTERISTIC SPEEDS 

A. Eigenvalue analysis 

A potential danger of zero-speed modes is that errors 
may couple to them, compound in place, and become 
inordinantly large. Zero-speed modes seem particularly 
hazardous in the case of moving black holes as they may 
adversely affect the dynamics by putting a drag on the 
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FIG. 3: Results of the gauge d t a = -2aK, dtfr = 
d t B l = t l - r]B i at times t = 4.5M (top panel) and t 



~i" > 

30M 



(bottom panel). The top panel clearly shows the lapse lagging 
behind the puncture, while the bottom panel exhibits noise 
propagating from the puncture. 



intractable. 

For spatial variation only in one dimension, with no 
initial transverse components, and assuming plane-wave 
solutions, the above system of equations can be written 
in the form: 



where 



d f \u) = ikU\ 



( a \ 

B 
b 

k 

h 
A 



(12) 



a i(kx— u)t) 



(13) 



with a, B, b, <f>, K, h, A, and T the constant amplitudes 
of a, B x , b x , <fi, K, h xx , A xx , and T x respectively, and 



puncture motion. In some cases zero-speed modes can 
also be a source of instability . So there is good reason 
to be aware of such modes. 

For the purpose of investigating characteristic speeds, 
a simple yet revealing method of analysis is to linearize 
the equations, assume plane wave solutions, and solve for 
the eigenmodes of the resulting algebraic system. With 
the gauge conditions described above, the BSSN equa- 
tions can be linearized about flat space as follows. 



(4) 
(5) 
(6) 



d = 


-2a K 


B l = 




V = 


\< B 


4> = 
k = 


~{aaK 

-didia + 


hij = 


-2a Aij 


Aij — 


—didja 



- diW) + p^ 
$d k K 



(7) 
(8) 



2a A tJ + d t V + 8,b l - -5 %] d k b k + ffidkhijQ) 



r = 



-2a didj 
4 
3 



-a d k dkhij + -a dif 3 

,iTF 



-aodjT 1 



f3 dkAij 
1 

3' 



(10) 

(11) 



where a = a — oiq, b l = fl 1 — /3q, hij = %j 



Sij, and 



the damping term on B % has been dropped as it should 
not affect the real part of the characteristic speeds. The 
lapse and shift are assumed to have constant, zeroth order 
terms ao and /3q, respectively which conveniently allows 
us to capture essential effects of the lapse factors and the 
advection terms without making the eigenvalue problem 
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The eigenval ues of M 
uo,-^Po + \\/Pl + 8a or 
\Po - 



are 0, 
\Po ~ 



-Po-Po + ap, 



-Po 
\Po 



4ag. 



8a 0> - 

If \v) is the 



yPl + 4ag, and -i/3 - \/pJ 
eigenvector associated with eigenvalue i>, and the covec- 
tor (v\ is defined such that J2 V \v)(v\ — I, then the eigen 
decomposition of M allows the system of evolution equa- 
tions to be written as a series of advection terms, each 
associated with a characteristic velocity equal to one of 
the eigenvalues: 



d t \u) = 

The inner product with < 

d t (v\u) 



E 



vd x \v){v\u) 



v\ gives: 

= v3t< 



(15) 



(16) 



Now we are interested in slow speed modes. Note 
that although P % is typically initialized at zero, with the 
Gamma-driver condition it evolves rapidly and signifi- 
cantly enough from zero that we do not classify Pq as a 
slow mode (since we're interested in moving punctures, 
P % is generally non- vanishing) . In addition to the obvious 
zero speed mode there are two modes that approach zero 
in the limit where ao approaches zero. If |ao| << \Po\ 
then -\Pq + \^Pl + %a and -\p + \\/ Pi + 4ag be- 
come approximately 4ao//3o an d 2«q/ Pq, respectively. In 
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the full nonlinear system, the 1+log slicing condition will 
collapse a nearly to zero in a finite region around the 
puncture (where /3 l will be nonzero for a moving punc- 
ture). Modes which propagate at a speed approximately 
proportional to a positive power of a for small a, which 
we will refer to as a-speed modes, are thus of potential 
concern. In particular, late in the evolution of black hole 
binaries there will emerge a significant, stationary region 
in which the lapse has collapsed nearly to zero, in which 
a-speed modes are effectively zero-speed modes. 

The eigencovector for v = — ^/3 + \y/ 0$ + 8ao in- 
volves only the fields a and K, is relatively independent 
from the other fields, and easy to fix, so we will address 
it first. This a-speed mode can be understood also in the 
context of the original nonlinear equations by noting that 
a couples with K to give a wave equation such that the 
speed of fluctuations of the lapse goes to zero when the 
lapse itself goes to zero. It is particularly egregious when 
the lapse couples to an a-speed mode, as once it collapses 
nearly to zero it will tend to get "stuck" there. But this 
can be remedied in a natural way simply by adding an 
advection term as follows, 

a = -2aK + (3 1 d j a, (17) 

which will modify the principal part of its wave equation 
such that when a goes to zero its speed will go to /3. Note 
that Eq. (|17|) . first used for moving punctures in Q, is 
consistent with the original Bona-Mass family of slicing 
conditions [ill]. 

The eigencovector for — i/3 + | v^o + ^ a o * s more 
complicated, involving not only a but several other fields 
as well. It turns out that advecting the lapse will not re- 
move this mode. But noting that this eigenfield involves 
b x , and that the equation for b x is conspicuously absent 
an advection term, an obvious stratagem is to advect the 
original j3 l as follows: 

P l = - A B l + pdjP l (18) 

Finally, Eq. I|16|) for the eigenvalue v = yields 

d t (B x -f x )=0 (19) 

This equation, a trivial consequence of the original 
Gamma-driver condition, immediately suggests a simple 
modification. The most natural way to remove this zero- 
speed mode is to advect it at (3 speed by modifying the 
B l equation thus: 

d t B l = d t r + ftdj (B l - P) - V B 1 (20) 

It turns out that the +p^djB z term alone will not suf- 
fice, as for example it might introduce an exponentially 
growing mode. Meanwhile the —j3 3 djT' t term alone adds 
another a-speed mode to the system. 

The foregoing analysis suggests the addition of four 
advection terms to the gauge equations, as indicated in 



Eqs. (|17I18I20|) . The lapse equation is relatively inde- 
pendent of the others: if it is advected then a particu- 
lar a-speed mode is removed, and if it is not advected 
then this particular a-speed mode remains. But the 
three advection terms suggested for the shift equations 
are strongly interdependent, so their combined effects are 
not so obvious and it is instructive to consider which un- 
desirable modes result from which combinations of these 
three terms. The possible combinations and their re- 
sulting "bad" speeds are summarized in Table 1. Here 
we see that for p > 0, only when all the above advec- 
tion terms are included will the system be free of slow- 
speed and exponentially growing modes; while for p = 0, 
the — ftdjT 1 term remains critical. In the case with all 
advection terms the benign eigenspeeds are: — (3q, —fio, 
-0 o -a o , -Po+a , -f3 -\/a%, -0o + V®o, -A— V 2 «o, 
— Po + \/2oq (where — (3q is listed twice because it is now 
associated with two distinct eigenvectors). 



Case # 


+P i d i B i 






"bad" speeds 


1 


N 


N 


N 


0, a p 


2 


N 


N 


Y 





3 


N 


Y 


N 


p/2 


4 


N 


Y 


Y 


a p 


5 


Y 


N 


N 





6 


Y 


N 


Y 


— i 


7 


Y 


Y 


N 




8 


Y 


Y 


Y 


none 



TABLE I: Effect of various advection terms in the shift equa- 
tions on the presence of undesirable eigenspeeds. Y' or 'N' 
indicates whether each advection term is included or not, re- 
spectively, in Eq. 11811 or Eq. 120H . 0, a p , or —i indicates 
whether the resulting linearized equations have a zero-speed 
mode, an a-speed mode, or an exponentially growing mode, 
respectively. The lapse is assumed to be advected as in 
Eq. 117H : otherwise an additional a-speed mode would ap- 
pear in every case. 



B. Numerical tests 

Evidence of the slow-speed modes found in the lin- 
earized analysis can often (but not always) be found in 
simulations of the full nonlinear system. Note that zero- 
speed modes only represent a potential danger. In some 
cases errors may not couple to these modes even though 
they exist in the equations. 

Here we provide a few examples of the cases given in 
Table U Fig. 0] represents Case #1 with p = 1, which 
is predicted to have a zero-speed mode and indeed man- 
ifests a non- propagating feature at x = — 3M. On the 
other hand, Case #1 with p = 0, also predicted to have 
a zero-speed mode, belies no indication of it in Fig. [SI 
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FIG. 4: Results for the gauge d t a = ~2aK + (3 j d 3 a, Q t fi = 
\aB\ d t B l = dtT* - r 1 B i (Case #1 in Table|T|with p=l) at 
time t = 14M. Non-propagating features are evident around 
x = — 3M. The conformal metric becomes singular by t = 
15M 

instead appearing quite smooth. In this case, evidently, 
the zero-speed mode does not couple significantly with 
the other fields. In Fig. El the spike in T x , which becomes 
steeper in time as the "tail" in T x grows, is apparently 
unrelated to slow-speed modes as none are predicted in 
this Case #3 with p = 0. However, for a longer evolution 
this growing gradient in T x can be expected to adversely 
affect the constraints. In Fig. the tail in T x and in 
particular the non-propagating bend in the tail around 
x = —2M seems to be evidence of the zero-speed mode 
predicted for Case #5. In Fig. [HI noise in T x grows ex- 
ponentially, as expected for Case #6. And finally, Case 
#8, depicted in Fig. ED and Fig. is demonstrated to 
be free of the previously identified "bad" speed modes 
whether p = or p = 1. (Although the former appears 
much smoother, perhaps because the eigenspeeds — /?o±l 
allow faster propagation of error away from the puncture 
than — (3q ± y/otQ.) 

IV. PROPERTIES OF THE "CLEANEST" 
GAUGES 

Of all the combinations suggested in Table [U two dis- 
tinguish themselves in numerical tests as conducive to 
particularly smooth propagation of the moving black 
hole. These gauges are Case #8 with p — 0, which we 
will refer to as the "shifting-shift case" : 

dt/P-ptdjF - ~B* (21) 

ihir I'iijr = d t r - d 3 r - V B l (22) 

and Case #1 with p — 0, which we will refer to as the 
"non-shifting-shift" case: 
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FIG. 5: Results for the gauge d t a = ~2aK + ftdja, dtfr = 
§B\ d t B* = d t r - TjB* (Case #1 in Table Q with p=0) at 
time t = 30M. The evolution appears smooth. 
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FIG. 6: Results for the gauge d t a = -2aK + P J d 3 a, dtft = 
\aB\ d t B l = <9 t r - ftdjf* - rjB* (Case #3 in TableQwith 
p=0) at time t = 30M. T x continues to grow larger and its 
gradient steeper. 

d t (3 l = - A B l (23) 

d t B l = dtF-riB* (24) 

The efficacy of Eqs. I|23I24|I in evolving moving punc- 
tures was originally demonstrated in |4J, Il7| . Apparently 
the zero-speed mode is not problematical in this case. 
Regarding the "shifting-shift" condition, Eqs. H21I22JI . we 
have recently used it successfully in the stable evolution 
of two equal mass black holes through 5.5 orbits. This 
latter gauge has also been recommended for the st rong 
hyperbolicity it brings to the BSSN equations 0, H9| . 
suggesting it may be a more robustly stable choice than 
Eqs. I|23I24|I . Nevertheless for the single black hole sim- 
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FIG. 7: Results for the gauge d t a = ~2aK + fidjOt, dtP* = 
\aB\ d t B i = dtT * + P j djB* - ryB* (Case #5 in TableQwith 
p=0) at time t = 30M. The feature around x = —2M does 
not propagate. 
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FIG. 8: Results for the gauge d t a = -2aK + P J d 3 a, dtp = 
\aB l + [3 j djf3\ d t B l = d t T l + /3 j d j B i - v B l (Case #6 in 
Table |3 with p=0) at time t = 10M The noise in f x grows 
exponentially. The conformal metric becomes singular by t = 
UM. 



ulations represented by Figs. ISland llUI the results from 
these two gauge choices are very similar. 

For these two gauge choices, we now consider varying 
the damping parameter, which has up to now been set 
to rj = 2. By moving toward r/ — > Eqs. (|21I22|) appear 
a little closer to realizing a "Gamma-freezing" condition. 
Indeed, with rj = 0, we find that the the black holes 
come closer to realizing the physically expected velocity. 
Both the "shifting-shift" and "non-shifting-shift" gauge 
options with rj = have proven to allow stable evolutions 
of a single black hole. In Fig. ^2 we find that in both cases 
setting rj — results in significantly smaller values for 
r s , meaning a closer approximation to the Dirac gauge. 
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FIG. 9: Results for the gauge d t a = ~2aK + ftdja, dtfr = 
^aB' + ^djF, d t B l = dtf l +f3 : >d 3 {B*-r)-r 1 B l (Case #8 
in TableQwith p=l) at time t — 30M. Aside from the sharp 
features in T l , which do not grow significantly in time, the 
propagation is non-pathological. 
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FIG. 10: Results for the gauge d t a = -2aK + /3 3 djCt, d t p i = 
|B" + P j d 3 f3\ d t B' = ftf' + /3 j dj(B i - f l ) - n B % (Case #8 
in TableQwith p=0) at time t = 30M. The evolution is very 
smooth. 



Further, we comment that we have observed numerically 
that for r) > 0, T l appears to show very slow linear growth 
in time, whereas for 77 = 0, T l appears to be bounded. 
Comparing the "shifting-shift" and "non-shifting-shift" 
options with 77 = we now find more noticeable differ- 
ences in r l , with our recommended "shifting-shift" op- 
tion giving a smoother result with T 1 holding closer to 
zero near the puncture. 

These two gauges share an additional feature in com- 
mon, which is that either Eqs. I23I24H or Eqs. (|21I22|) 
can be integrated to give the relation 

B l = f 1 - ^ (25) 
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FIG. 11: F x for the gauge d t a = -2aK + ft dja, d t f3 i = §B\ 
9t-B* = dtT 1 (Case #1 with rj — 0) and for the gauge fta = 
-ZaK+FdjOt, dtP* = f B'+FdjP 1 , d t B l = dtr+ftdjiB 1 - 
Y i ) (Case #8 with 77 = 0) at time t = 10M. In the "shifting- 
shift" case (Case #8) F 21 is smaller and smoother around the 
puncture than in the "non-shifting-shift" case (Case #1). 
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FIG. 13: Results for the gauge d t a = -2aK + ftdja, dtp = 
§f i + ftdjft - rjff at time t = 30M. The results are nearly 
identical to those depicted in Fig. I1UI 



term. Dispensing with the B l evolution equation results 
in a slightly more efficient numerical implementation, and 
also guarantees removal of the zero-speed mode associ- 
ated with B l . 

Finally it may be noted that when rj = and r* — » 0, 
Eq. HIJ) (and Eqs. 112112211 ) admits a "shock" solution of 
the form (x — xo)/(to — t). This represents the nonlinear 
tendency of this equation to advect large magnitudes of 
the shift faster than smaller magnitudes, which for neg- 
ative values and a positive slope can lead to a vertical 
slope. However, as dtT 1 depends on derivatives of /?% T z 
is unlikely to vanish when those derivatives become large. 
Numerically (3 l has proven to be very well behaved so T l 
evidently acts as an effective "shock absorber", at the 
expense of not vanishing entirely. 

V. CONCLUSIONS 



FIG. 12: Results for B x and f x - §77/3* with the gauge d t a = 
-2aK+Fd ja , dtP 1 = \B l + Fd 3 P\ d t B l = dtP+ftdjiB'- 
T ! ) — r/B 1 at time t = IBM. The difference between the two 
curves is within 10 -11 . 

since B l — T l = j3 l = initially. This fact, also evi- 
dent numerically (Fig. IT2T). suggests substituting for B % 
in the evolution equation for (3 l , to obtain (in the case of 
Eqs. (ED) 

dt(3 l = ^P+pdjF-riF (26) 

Fig. El demonstrates that the resulting shift condition, 
Eq. H26fl . can be evolved numerically to yield the same 
stable, smooth simulation as with the analytically equiv- 
alent condition, Eqs. I|21I22|I . used previously. Similar 
success should also be obtainable without the ftdjft* 



We have investigated gauge conditions that are appro- 
priate for the moving puncture approach to black hole 
simulations. Potential hazards from zero-speed or slow- 
speed modes have been identified, a methodology for an- 
alytically exploring various gauge choices has been pre- 
sented, and a gauge free of slow-speed modes has been 
recommended. Several gauges were also studied numer- 
ically and two gauges distinguished themselves as par- 
ticularly well-adapted to smooth black hole motion. In 
both cases we have suggested a simplification of the shift 
evolution equations, as well as explored the possibility of 
eliminating a traditional damping term to better realize a 
"Gamma-freezing" condition. We found that eliminating 
this damping term on the shift does indeed minimize T l 
and yield smooth evolutions, more so with the addition 
of our recommended "shifting-shift" terms, and warrants 
further experimentation. We intend to investigate this 
rj — > limit more thoroughly in future work. 



8 



Acknowledgments 

This work was supported in part by NASA grants 
ATP02-0043-0056 and O5-BEFS-05-0044. The simula- 
tions were carried out using Project Columbia at NASA 
Ames Research Center and at the NASA Center for 
Computational Sciences at Goddard Space Flight Center. 



M.K and J.v.M. were supported by the Research Associ- 
ateship Programs Office of the National Research Coun- 
cil and the NASA Postdoctoral Program at the Oakridge 
Associated Universities. We are very grateful to Carsten 
Gundlach for carefully reading and pointing out an error 
in an earlier draft. 



[1] B. Bruegmann, W. Tichy, and N. Jansen, Phys. Rev. 

Lett. 92, 211101 (2004), gr-qc/0312112. 
[2] M. Alcubierre et al., Phys. Rev. D72, 044004 (2005), 

gr-qc/0411149. 

[3] J. G. Baker, J. Centrella, D.-I. Choi, M. Koppitz, and 
J. van Meter, Phys. Rev. Lett. 96, 111102 (2006), gr- 
qc/0511103. 

[4] M. Campanclli, C. O. Lousto, P. Marronetti, and 
Y. Zlochower, Phys. Rev. Lett. 96, 111101 (2006), gr- 
qc/05 11048. 

[5] J. G. Baker, J. Centrella, D.-I. Choi, M. Koppitz, and 
J. van Meter, Phys. Rev. D73, 104002 (2006), gr- 
qc/0602026. 

[6] T. Nakamura, K. Oohara, and Y. Kojima, Prog. Theor. 

Phys. Suppl. 90, 1 (1987). 
[7] M. Shibata and T. Nakamura, Phys. Rev. D 52, 5428 

(1995). 

[8] T. W. Baumgarte and S. L. Shapiro, Phys. Rev. D 59, 

024007 (1999), gr-qc/9810065. 
[9] B. Imbiriba et al., Phys. Rev. D70, 124025 (2004), gr- 

qc/0403048. 



[10] M. Alcubierre, B. Briigmann, P. Diener, M. Koppitz, 

D. Pollney, E. Seidel, and R. Takahashi, Phys. Rev. D 

67, 084023 (2003), gr-qc/0206072. 
[11] C. Bona, J. Masso, E. Seidel, and J. Stela, Phys. Rev. 

Lett. 75, 600 (1995), gr-qc/9412071. 
[12] S. Brandt and B. Briigmann, Phys. Rev. Lett. 78, 3606 

(1997), gr-qc/9703066. 
[13] J. D. Brown and L. L. Lowe, J. Comput. Phys. 209, 582 

(2005), gr-qc/0411112. 
[14] P. A. M. Dirac, Phys. Rev. 114, 924 (1959). 
[15] L. Smarr and J. W. York, Phys. Rev. D17, 1945 (1978). 
[16] M. Alcubierre, G. Allen, B. Briigmann, E. Seidel, and 

W.-M. Suen, Phys. Rev. D 62, 124011 (2000), gr- 

qc/9908079. 

[17] M. Campanelli, C. O. Lousto, and Y. Zlochower, Phys. 

Rev. D 73, 061501 (2006), gr-qc/0601091. 
[18] C. Gundlach and J. M. Martin-Garcia (2006), gr- 

qc/0604035. 

[19] H. Beyer and O. Sarbach, Phys. Rev. D70, 104004 
(2004), gr-qc/0406003. 



